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Accuracte calculations on the bond length, dipole moment and harmonic 
frequency of CO are presented, using large basis sets and high levels of electron 
correlation. The geometric structure, forceconstants and binding energies of 
Cr(C0)6 and Cr(C0)5 are computed using large basis sets and high levels of 
electron correlation. The molecule Oj is studied, using large basis sets and high 
levels of electron correlation, including CASSCF, CASSI and CASPT2 methods. 
Binding energies, geometries and frequencies are computed. Symmetry breaking is a 
particular problem for the antisymmetric stretch, which is addressed using the 
CASSI method. The symmetry breaking problem in O4 has also been studied using 
the Brueckner coupled-cluster method. This gives results in good agreement with 
CASSI. A multi-region numerical integration scheme is investigated for use in 
Density Functional Calculations. This scheme is found to give comparable results to 
a widely used scheme based on the Euler- Maclaurin technique. 



“Bond length, dipole moment and harmonic frequency of CO”, L. A. Barnes, 
B. Liu and R. Lindh, J. Chem. Phys., 98, 3972-3977, (1993) 

A detailed comparison of some properties of CO is given, at the modified 
coupled-pair functional (MCPF), single and double excitation coupled-cluster 
(CCSD) and CCSD(T) levels of theory (including a perturbational estimate 
for connected triple excitations), using a variety of basis sets. With very large 
one-particle basis sets, the CCSD(T) method gives excellent results for the 
bond distance, dipole moment and harmonic frequency of CO. In a 
[6s 5p 4d 3/ 2 g l/i] + (Is lp Id) basis set, the bond distance is about 0.005 ao 
too large, the dipole moment about 0.005 a.u. too small and the frequency 
about 6 cm -1 too small, when compared with experimental results. 

“Structure and energetics of Cr(CO)6 and Cr(CO)s”, L. A. Barnes, B. Liu and 
R. Lindh, J. Chem. Phys., 98, 3978-3989, (1993) 

The geometric structure of Cr(CO)6 is optimized at the modified coupled-pair 
functional (MCPF), single and double excitation coupled-cluster (CCSD) and 
CCSD(T) levels of theory (including a perturbational estimate for connected 
triple excitations), and the force constants for the totally symmetric 
representation are determined. The geometry of Cr(CO)s is partially 
optimized at the MCPF, CCSD and CCSD(T) levels of theory. Comparison 
with experimental data shows that the CCSD(T) method gives the best 



results for the structures and force constants, and that remaining errors are 
probably due to deficiencies in the one-particle basis sets used for CO. A 
detailed comparison of the properties of free CO is therefore given, at both the 
MCPF and CCSD/CCSD(T) levels of treatment, using a variety of basis sets. 
With very large one-particle basis sets, the CCSD(T) method gives excellent 
results for the bond distance, dipole moment and harmonic frequency of free 
CO. The total binding energies of Cr(C0)6 and Cr(C0)5 are also determined 
at the MCPF, CCSD and CCSD(T) levels of theory. The CCSD(T) method 
gives a much larger total binding energy than either the MCPF or CCSD 
methods. An analysis of the basis set superposition error (BSSE) at the 
MCPF level of treatment points out limitations in the one-particle basis used 
here and in a previous study. Calculations using larger basis sets reduce the 
BSSE, but the total binding energy of Cr(CO)e is still significantly smaller 
than the experimental value, although the first CO bond dissociation energy 
of Cr(C0)6 is well described. An investigation of 3s3p correlation reveals only 
a small effect. The remaining discrepancy between the experimental and 
theoretical total binding energy of Cr(C0)6 is probably due to limitations in 
the one-particle basis, rather than limitations in the correlation treatment. In 
particular an additional d function and an / function on each C and 0 are 
needed to obtain quantitative results. This is underscored by the fact that 



even using a very large primitive set (1042 primitive functions contracted to 
300 basis functions), the superposition error for the total binding energy of 
Cr(CO)6 is 22 kcal/mol at the MCPF level of treatment. 

• “The fraternal twins of quartet Oj”, R. Lindh and L. A. Barnes, J. Chem. 
Phys., 100, 224-237, (1994) 

Eleven stationary geometries of quartet Oj have been studied by ab initio 
methods. The geometries were optimized at the CASSCF level of theory and 
the energies were calculated by the CASPT2 method, using DZP, TZ2P, 
ANO[5s4p2d] and ANO [6s5p3d2/] basis sets. The rectangular and 
trans-planar structures are found to be the most stable, with an energy barrier 
to conversion between the two at the threshold of dissociation. Both have a 
delocalized hole and are stable relative to separated O2 and 0 J by 11.0 and 
11.5 kcal/mol for the rectangular and the trans-planar structure, respectively, 
compared with the experimentally deduced energy in the range of 9.2 to 10.8 
kcal/mol. The adiabatic ionization potentials of O4 and O2 are computed to 
be 11.67 and 12.21 eV, while experimental values are 11.66 and 12.07 eV, 
respectively. The vibrational frequencies have been computed for all degrees of 
freedom at the CASSCF level of theory. Symmetry breaking is found to be a 
particular problem in the computation of the antisymmetric stretch frequency 
for the delocalized structures at the CASSCF level of theory. Attempts to 



rectify these problems using the RASSCF method leads to additional 
difficulties, but further analysis yields insight into the symmetry breaking and 
problems with earlier calculations. Finally, a non-orthogonal Cl calculation 
based on the interaction of localized CASSCF wavefunctions using the CASSI 
method leads to a balanced treatment of the antisymmetric stretch which is 
free from symmetry breaking. The study explains the four most prominent 
absorption frequencies observed in the partially unassigned IR spectrum of O 4 
isolated in solid neon as the antisymmetric OO-stretch, and the combination 
band of the symmetric and antisymmetric OO-stretch of both the rectangular 
and frans-planar structures. 

• “Symmetry Breaking in 0| : an application of the Brueckner Coupled Cluster 
Method”, L. A. Barnes and R. Lindh, Chem. Phys. Lett., 223 , 207, (1994) 

A recent calculation of the antisymmetric stretch frequency for the rectangular 
structure of quartet O 4 using the QCISD(T) method gave a value of 
3710 cm -1 . This anomalous frequency is shown to be a consequence of 
symmetry breaking effects, which occur even though the QCISD(T) solution 
derived from a delocalized SCF reference function lies energetically well below 
the two localized (symmetry-broken) solutions at the equilibrium geometry. 
The symmetry breaking is almost eliminated at the CCSD level of theory, but 
the small remaining symmetry breaking effects are magnified at the CCSD(T) 



level of theory so that the antisymmetric stretch frequency is still significantly 
in error. The use Brueckner coupled cluster method, however, leads to a 
symmetrical solution which is free of symmetry breaking effects, with an 
antisymmetric stretch frequency of 1322 cm -1 , in good agreement with our 
earlier calculations using the CASSCF/CASSI method. 

• “A multi-region integration scheme”, L. A. Barnes, work in progress (to be 
published) 

In this preliminary report, a multi-region radial integration is compared to the 
recently proposed method due to Handy et al.. Preliminary results for small 
systems indicate that the new integration scheme is generally comparable to 
and sometimes better than that of Handy et al., although this conclusion is by 
no means firm. Work for larger systems is continuing. 



Symmetry breaking in O 4" : an application of the 
Brueckner coupled- cluster method 
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Abstract 

A recent calculation of the antisymmetric stretch frequency for the rectangular structure 
of quartet 0| using the QCISD(T) method gave a value of 3710 cm" 1 . This anomalous 
frequency is shown to be a consequence of symmetry breaking effects, which occur even 
though the QCISD(T) solution derived from a delocalized SCF reference function lies 
energetically well below the two localized (symmetry-broken) solutions at the equilibrium 
geometry. The symmetry breaking is almost eliminated at the CCSD level of theory, but 
the small remaining symmetry breaking effects are magnified at the CCSD(T) level of 
theory so that the antisymmetric stretch frequency is still significantly in error. The use 
Brueckner coupled cluster method, however, leads to a symmetrical solution which is free 
of symmetry breaking effects, with an antisymmetric stretch frequency of 1322 cm , in 
good agreement with our earlier calculations using the CASSCF/CASSI method. 

"Mailing address: NASA Ames Research Center, Moffett Field, California 94035-1000 
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Introduction 


Recently, Peel [1] reported ab initio calculations on some high symmetry isomers of quar- 
tet O4 using the singles and doubles quadratic configuration interaction method with a 
perturbational estimate of connected triple excitations (QCISD(T)) [2j. For the rectan- 
gular structure, he obtained an antisymmetric stretch frequency of 3710 cm -1 , which he 
noted was “unphysical”, but nevertheless was used to give a large zero-point correction 
to the energy of this isomer. 

Previously [3], we had carried out a detailed study of several structures of quartet O4 , com- 
puting amongst other things vibrational frequencies, relative energies and isotopic shifts. 
The methods used included complete active space self consistent field (CASSCF) [4], 
restricted active space self consistent field (RASSCF) [5], a second order perturbation 
method based on a CASSCF reference function (CASPT2) [6] and a non-orthogonal Cl 
method based on non-orthogonal CASSCF or RASSCF solutions (the complete/restricted 
active space state interaction (CASSI/RASSI) method [7]), using large generally con- 
tracted atomic natural orbital basis sets. These methods gave good results and prompted 
a new analysis [8] of the experimental vibrational spectrum of Oj , which supported our 
assignment. Of particular importance for the current work is our detailed analysis of 
symmetry breaking effects in the calculation of the antisymmetric stretch frequency for 
quartet 0|. We showed that spurious frequencies for quartet Oj were due to symme- 
try breaking effects due to the competition between localized and delocalized structures 
— we refer the reader to our original work [3], where we also give many references to 
earlier work on symmetry breaking. In our case, the symmetry breaking effects were re- 
solved through the use of the CASSI method, allowing the two non-orthogonal “localized” 
CASSCF wavefunctions to interact and give the correct qualitative form for the potential 
energy surface. 

On the basis of the form of the potential curve for the antisymmetric stretch around 
the equilibrium point for the rectangular structure, Peel [1] concluded that no symme- 
try br eak in g effects were evident at the QCISD(T) level of theory. In the first part of 
this work we demonstrate that this is incorrect: symmetry breaking effects are entirely 
responsible for the anomalous antisymmetric stretch frequency of 3710 cm -1 . For con- 
sistency, we use the same methods as in Ref. [1] — QCISD(T) calculations based on an 
unrestricted Hartree Fock (UHF) reference function and a 6-31G* basis set. The geometry 
(Roo=1.186 A, Rcm=2.378 A) was also taken from Ref. [1]. Roo is the intra-fragment 
0-0 bond distance, Rcm is the inter-fragment bond distance (the distance between the 
center of masses of the two fragments). The rectangular structure is illustrated in Fig. 1(a), 
which also shows Roo and Rcm- 

Following this, we give the results of calculations at the CCSD(T) level of theory (singles 
and doubles coupled-cluster plus a perturbational estimate of the effects of connected 
triples excitations [9]). The geometry is reoptimized and from the computed potential 
curves and frequencies it is evident that symmetry breaking effects also occur at the 
CCSD(T) level of theory. Brueckner coupled cluster theory [10]-[15] has previously been 
used to treat symmetry breaking effects at the coupled cluster level of theory [16]. We 
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have applied this approach to the rectangular structure of quartet O 4 , using Brueckner 
coupled cluster calculations which include a perturbational estimate of the effects of con- 
nected triples excitations (BD(T)) [14], with semi-canonical orbitals [17]. Since we are 
interested in the qualitative nature of the results, rather than strictly quantitative results, 
we continue to use the 6-3 1 G* basis set and UHF reference function, rather than the larger 
generally contracted basis sets of our earlier study [3]. We note that in the CCSD(T) and 
BD(T) calculations, all electrons are correlated, whereas in the QCISD and QCISD(T) 
calculations, the Is core electrons were not correlated. 

The calculations were performed with the ACES II 1 suite of programs using IBM RISC 
SYSTEM/6000 computers at NASA Ames Research Center. 


Results and Discussion 

As discussed in detail in our previous work [3], the symmetry breaking in 0| is mani- 
fested by the existence of a delocalized solution, which exhibits D 2 h symmetry, and two 
localized solutions exhibiting C 2v symmetry which are mirror images of each other, having 
equal energies at the symmetric (D 2 *) point. The structures are illustrated in Fig. 1, with 
Fig. 1(a) showing the symmetric geometry, and Figs. 1 (b) and (c) showing the two symme- 
try broken localized solutions. In each case, the shorter 0-0 bond distance corresponds 
closely to the bond distance of Oj, whereas the longer bond distance corresponds closely 
to the bond distance of 0 2 [3]. Thus the positive charge is localized on the bottom 0 2 
unit in Fig. 1 (b), and on the top 0 2 unit in Fig. 1 (c). The symmetry breaking vibrational 
mode is the intra-fragment antisymmetric stretch u> 5 , and this is illustrated in Figs. 1 (b) 
and (c) by the arrows. 

In Table I, we give the total energies of the different solutions at the UHF , QCISD and 
QCISD(T) levels of theory. Comparing the total energy at the QCISD(T) level with 
that given in Ref. [1], we see that the results of Ref. [1] are based on the delocalized 
reference function. It is interesting to compare the energy differences of the delocalized 
and localized solutions at the various levels of theory. At the UHF and QCISD levels 
of theory, the localized solution is below the delocalized solution by about 3 kcal/mol. 
However, at the QCISD(T) level of theory, the localized solution is more than 7 kcal/mol 
higher in energy than the delocalized solution. Thus it seems that the D 2 h structure is 
favoured at the highest level of theory. 

To understand the origin of the spurious frequency at the QCISD(T) level of theory, we 
have computed the energy as a function of the antisymmetric stretch coordinate (see 

iACES II is a computational chemistry package especially designed for coupled cluster and many body 
perturbation calculations. The SCF, transformation, correlation energy and gradient codes were written 
by J. F. Stanton, J. Gauss, J. D. Watts, W. J. Lauderdale and R. J. Bartlett. The two-electron integrals 
are taken from the vectorized MOLECULE code of J. Almlof and P. R. Taylor. ACES II includes a 
modified version of the ABACUS integral derivatives program, written by T. Helgaker, H. J. Jensen, P. 
J0rgensen, J. Olsen, and P. R. Taylor, and the geometry optimization and vibrational analysis package 
written by J. F. Stanton and D. E. Bernholdt. 
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also Refs. [1] and [3]). The results are presented in Figs. 2-4, where we have used the 
same scales in order to facilitate comparison of the different curvatures resulting from 
the different methods. The delocalized solution is given by Ed e i and the two localized 
solutions are given by E a and Eb in each case. Also given in Figs. 3 and 4 are curves at 
the CCSD, CCSD(T), BD and BD(T) level of theory, which we discuss later. Using the 
notation of Figs. 1(b) and (c), AR is defined as Ri— R 2 . 

The behaviour of the energy at the correlated level of theory is driven by the behaviour 
of the UHF energies, given in Fig. 2. The delocalized solution lies above the localized 
solutions, and as we follow the antisymmetric coordinate, the delocalized solution rapidly 
approaches the localized solutions, until at ARwO.0185 A the localized and delocalized 
solutions merge. Thus the QCISD and QCISD(T) energies are constrained by this fact 
— the delocalized energy and one of the localized energies (depending on whether the 
distortion is positive or negative) must be equal at ARft*0.0185 A. In Fig. 3, we see that 
at the QCISD level of theory the energies are significantly better in the sense that the 
localized solution energies are much flatter (tending towards delocalized solutions). In 
this case the flatness of the localized solutions means that the delocalized solution is 
driven down in energy in order to meet the constraint of equal energy at AR~0.0185 A, 
leading to an imaginary frequency at the QCISD level of theory. At the QCISD(T) level of 
theory (Fig. 4) the delocalized solution is well below the localized solution, and so is driven 
rapidly upward to meet the equal energy constraint, resulting in the very large 3710 cm -1 
frequency. Thus the QCISD method does not entirely overcome the inherent problems 
with the UHF reference function, and the triples perturbation correction is unable to 
overcome the residual problems with the QCISD method. We note that we found similar 
problems with the CASPT2 method, which was not able to overcome the problems of a 
localized CASSCF reference function [3]. 

In Tables II and III we present the results from the CCSD(T) calculations. From the 
energy separations, we see that while the UHF separation is very similar to that given in 
Table I (which has a slightly different geometry), the CCSD and CCSD(T) separations 
are very different to those at the QCISD and QCISD(T) levels of theory. The delocal- 
ized and localized solutions are very close in energy at the CCSD level of theory, and 
unlike the QCISD results the localized solution is above the delocalized solution. As for 
the QCISD(T) results, the perturbational triples correction increases the separation be- 
tween the delocalized and localized solutions at the CCSD(T) theory when compared with 
CCSD, although the effect is much smaller than the QCISD and QCISD(T) difference. 
We note that one difference between the QCISD and CCSD calculations was that the Is 
core electrons were excluded from the calculations at the QCISD and QCISD(T) levels of 
theory, whereas they were included in the calculations at the CCSD and CCSD(T) levels 
of theory. To check whether this difference has any effect on the symmetry breaking at the 
QCISD/QCISD(T) level of theory, we also computed the separation between the delocal- 
ized and localized solutions at the symmetric point including the Is electrons using these 
methods. The separations are barely different from the original results, so we conclude 
that removing the core Is electrons from the QCISD calculations is not the cause of the 
large difference between the CCSD/CCSD(T) and QCISD/QCISD(T) results. 

In Figs. 3 and 4 we present the CCSD and CCSD(T) potential curves for the antisymmetric 
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stretch, which may be compared with the QCISD and QCISD(T) curves on the same 
figures. The behaviour of the UHF reference function energies around the CCSD(T) 
equilibrium geometry for the antisymmetric stretch is very similar to that given in Fig. 2 
around the QCISD(T) equilibrium geometry, so we may discuss the CCSD and CCSD(T) 
curves in the same light as the QCISD and QCISD(T) curves. Thus the CCSD and 
CCSD(T) energies have similar constraints to the QCISD and QCISD(T) energies — 
irrespective of the separation at AR=0 (the symmetric geometry), at AR«0.0185A the 
delocalized and one of the localized energies must be equal. Inspection of the curves shows 
that this is so. However, the most striking difference between the CCSD and QCISD 
curves comes from the fact that the curves at the CCSD level are much closer together, 
so that this constraint has only a small effect on the antisymmetric stretch frequency at 
the CCSD level of theory. In fact, the antisymmetric stretch frequency at the CCSD level 
of theory is a very reasonable 1220 cm" 1 (at the CCSD equilibrium geometry), which is 
to be compared with a value of around 1500? cm' 1 at the QCISD level of theory. Overall, 
we see that the CCSD approach has almost eliminated the symmetry breaking effects. 

As discussed above, the addition of the perturbative triples correction increases the sep- 
aration between the delocalized and localized solutions at the symmetric point, and this 
is evident in Fig. 4. Thus the delocalized curve at the CCSD(T) level is more affected 
by symmetry breaking than the CCSD curve, although this effect is much smaller than 
that found with the QCISD(T) method. Thus the antisymmetric stretch frequency at 
the CCSD(T) level is 1922 cm -1 , compared with 3710 cm -1 at the QCISD(T) level of 
theory. The origin of this difference is quite evident from the potential curves — it is the 
large difference in separations at the symmetric point. The other remarkable feature of 
the CCSD(T) potential curves is the near coincidence of the two localized curves, which 
is again quite different to the QCISD(T) results. Thus the CCSD(T) approach is quite 
close to removing the symmetry breaking effects, but is still not able to overcome the 
small deficiencies evident at the CCSD level of theory. 

The geometry at the CCSD(T) level (Table III) is very similar to that found at the 
QCISD(T) level of theory [1], and for the most part the frequencies are quite similar to 
those given in Ref. [1]. The exceptions are the antisymmetric stretch u> 5 (discussed above) 
which changes from 3710 to 1922 cm" 1 , and the (inter-fragment) antisymmetric stretch 
u> 6 which is reduced from 595 to 97 cm" 1 . The CCSD(T) value for a* is in accord with our 
earlier results [3] and the results for the frans-planar structure [1,3]. Thus it seems that 
the QCISD(T) value for u 6 is significantly too high also. Considering Fig. 1, the mode u> 6 
may be envisioned in an analogous way to u>5, except that the distortion occurs along the 
Rcm direction instead of the Roo direction. Thus it is possible, though less likely (due 
to the large inter-fragment bond distance Rcm)> f° r symmetry breaking effects to occur 
for u>6 also. This would involve localization on the left and right sides of O4 rather than 
the top and bottom, which occurs for u; 5 . However, we have not investigated this in any 

detail here. 

The results at the BD(T) level of theory axe given in Table IV. The geometry optimized 
at the BD(T) level of theory is the same as that of the CCSD(T) level of theory (which 
was constrained to have D 2 h symmetry, whereas the BD(T) calculation was not), and 
the BD(T) energy is also very similar to the CCSD(T) energy. To investigate whether 
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symmetry breaking effects are still present at the BD(T) level of theory, we have again 
plotted the energy as a function of the antisymmetric distortion, and the results are given 
in Figs. 3-5. 

As we saw previously, the behaviour of the correlated methods was constrained by the 
behaviour of the reference function. In Fig. 5 we give the Brueckner reference determi- 
nant energy and the UHF energy from which the Brueckner calculation was initiated. 
Before discussing these results, we emphasize that the comparison between the Brueckner 
reference determinant energy and the UHF energy is not rigourous since the Brueck- 
ner reference determinant is a product of the correlated calculation. Nevertheless, it is 
enlightening. 

The UHF solutions of Fig. 5 are (qualitatively speaking) a subset of those given in Fig. 2. 
The curve is discontinuous because we varied AR with a larger stepsize than for Fig. 2, 
and the SCF converged to solutions on different potential curves at different points, rather 
than the solution on the same potential curve as in Fig. 2. The character of the UHF 
orbitals is of course very different for the different potential energy surfaces, varying from 
delocalized to localized on the top of the molecule or localized on the bottom of the 
molecule, and this variation is reflected in the energies. In contrast to this, the Brueckner 
reference energy is very smooth despite the large changes in the UHF orbitals from which 
it began, indicating that the Brueckner approach is not affected by the starting orbitals. 
At the symmetric point we have also verified that the Brueckner method is independent 
of the starting orbitals — whether localized or delocalized UHF orbitals are used, the 
Brueckner approach leads to the same symmetric (delocalized) solution. Thus there is only 
one solution at the Brueckner level of theory. This behaviour is in accord with previous 
studies [16] using the Brueckner approach for other systems which exhibit symmetry 
breaking. 

At the BD and BD(T) levels of theory (Figs. 3 and 4) we see that the antisymmetric 
stretch is very smooth and gives a positive frequency, which is 1322 cm 1 at the BD(T) 
level of theory (Table IV). It is interesting to compare the different curvatures for the 
different methods in Figs. 3 and 4. It is evident that the CCSD and CCSD(T) curvatures 
are much closer to the BD and BD(T) curves than are those from QCISD and QCISD(T). 

We note from our previous study [3] that there is a significant basis set effect for the 
antisymmetric stretch frequency. At the CASSI level using a TZ2P basis set, the an- 
tisymmetric stretch was 1271 cm -1 whereas an ANO[5s4p2d] basis set gave a value of 
1259 cm" 1 and an ANO(6s5p3d2/] basis gave a value of 1296 cm' 1 . Considering the fact 
that the BD(T) approach should give a larger proportion of the dynamical correlation 
energy than our earlier frequency calculations at the CASSCF/CASSI level of theory, the 
agreement between the BD(T) frequency and our earlier values is very good. Thus the 
BD(T) results are very encouraging and in a large one particle basis this method should 
give very accurate results. In our earlier work [3] we showed that the dipole derivative at 
the CASSI level of theory was unphysically high. It would be of some interest to compute 
this quantity at the BD(T) level (in a large one particle basis) to determine whether a 
more reasonable dipole derivative would be obtained. 

In Table IV we also give the symmetric stretch frequencies at the BD(T) level of theory. 


6 



Given the agreement between the geometries at the CCSD(T) and BD(T) levels of theory, 
it is not surprising that the symmetric stretch frequencies are very similar for the two 
methods (and also in good agreement with the CASSCF results [3]). These results also 
support our earlier isotopic substitution analysis [3], where we used the CASSCF frequency 
for the symmetric stretch and the CASSI frequency for the antisymmetric stretch. 

To conclude, the antisymmetric stretch of quartet O 4 is significantly affected by symmetry 
breaking. As we discussed previously [3], it is necessary to properly account for this before 
a reliable frequency can be obtained. In the current work we have shown in detail how the 
previous [1] antisymmetric stretch frequency at the QCISD(T) level of theory is affected 
by symmetry breaking so that any analysis of the relative energies of the rectangular 
and trans-planar structures which includes zero-point corrections based on this frequency 
must be significantly in error. The CCSD approach gives significantly better results 
than QCISD, almost eliminating the symmetry breaking effects. However, the small 
remaining symmetry breaking effects are magnified at the CCSD(T) level of theory, so 
that the antisymmetric stretch is still affected significantly at the CCSD(T) level. The 
Brueckner coupled- cluster method (BD(T)), however, eliminates the symmetry breaking 
effects entirely, giving a single symmetric solution with an antisymmetric stretch frequency 
in good agreement with our earlier result at the CASSCF/CASSI level of theory [3]. This 
must make the BD(T) approach the method of choice for very accurate calculations when 
symmetry breaking is a potential problem and more than just a few electrons must be 
correlated. 
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Table I: Total energies (in Hartree) and energy differences (in kcal/mol) of quartet O 4 at 
the rectangular (D 2 /i) geometry, for the delocalized and localized solutions. The geometry 
is taken from the QCISD(T) calculations of Ref.[l] (Roo = l-186 A,Rcm=2.378 A) 


Method 

Edel 

Eloc 

AE 

UHF 

-298.7398444 

-298.7439668 

2.59 

QCISD 

-299.4774306 

-299.4830304 

3.51 

QCISD(T) 

-299.5050169 

-299.4933392 

-7.33 


Table II: Total energies (in Hartree) and energy differences (in kcal/mol) of quartet O 4 at 
the rectangular (D 2/l ) geometry, for the delocalized and localized solutions. The geometry 
is from the CCSD(T) approach, given in Table III 


Method 

Edel 

Eloc 

AE 

UHF 

-298.7406630 

-298.7447092 

2.54 

CCSD 

-299.4838322 

-299.4832236 

-0.38 

CCSD(T) 

-299.5127862 

-299.5107674 

-1.27 
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Table III: Total energy (in Hartree), geometrical parameters (in A) and frequencies (in 
cm -1 ) of quartet 0| at the rectangular (D 2 /i) geometry, using the CCSD(T) approach. 


E 

-299.5127862 

Roo 

1.1846 

Rcm 

2.3751 

^i(a s ) 

1713 

u; 2 (a s ) 

271 

w 3(bi s ) 

372 

<*>4(a„) 

175 

W5(b 3u ) 

1922° 

k>6(b 2u ) 

97 


0 This frequency was computed as a finite difference of energies, since the symmetry 
breaking effects lead to spurious results from the automatic finite difference of gradients 
approach of Aces II (which used different UHF solutions at different points of the finite dif- 
ference procedure). However, the automated approach was used for the other frequencies, 
including u)e (which was also computed by finite difference of energies as a check). 


Table IV: Total energy (in Hartree), geometrical parameters (in A) and frequencies (in 
cm -1 ) of quartet 0| at the rectangular (D 2 h) geometry, using the BD(T) approach. 


E 

-299.5127780 

Roo 

1.1846 

Rcm 

2.3751 

wi(a») 

1713 

w 2 (a fl ) 

270 

^ 5 ( 63 ,,) 

1322 
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Figure captions 


Figure 1 : Schematic structures of rectangular 0| . (a) The symmetric structure (D 2 * 
symmetry), (b) The bottom localized structure (C 2 „ symmetry), (c) The top localized 
structure (C 2v symmetry). The antisymmetric stretch mode tu 5 is illustrated 
schematically with arrows. 

Figure 2: UHF energies (around the QCISD(T) equilibrium geometry) as a function of 
the antisymmetric stretch coordinate for the rectangular structure of O 4 ". E a is given by 
O, Eb is given by + and Ed e i is given by □. 

Figure 3: QCISD, CCSD and BD energies (around the QCISD(T), CCSD(T) and 
BD(T) equilibrium geometries, respectively) as a function of the antisymmetric stretch 
coordinate for the rectangular structure of O 4 . The QCISD energies E a are given by O, 
Eb are given by + and Ea e i are given by □. The analogous CCSD energies are given by 
x, A and *, respectively. The BD energies are given by,. For clarity, the CCSD energies 
are offset by -0.005 Hartree, whereas the BD energies are offset by -0.01 Hartree. 

Figure 4: QCISD(T), CCSD(T) and BD(T) energies as a function of the antisymmetric 
stretch coordinate for the rectangular structure of O 4 . The QCISD(T) energies E a are 
given by O, Eb are given by -f and Edei are given by □. The analogous CCSD(T) 
energies are given by x, A and *, respectively. The BD(T) energies are given by,. For 
clarity, the CCSD(T) energies are offset by 0.0025 Hartree. 

Figure 5: UHF and Brueckner determinant reference energies as a function of the 
antisymmetric stretch coordinate for the rectangular structure of O 4 . The UHF energies 
are given by O and the Brueckner determinant reference energy is given by + 
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Abstract 

In this preliminary report, a multi-region radial integration is compared to 
the recently proposed method due to Handy et a/.. Preliminary results for small 
systems indicate that the new integration scheme is generally comparable to 
and sometimes better than that of Handy et a/., although this conclusion is by 
no means firm. Work for larger systems is continuing. 
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1 Introduction 


There is recently a great deal of interest in the testing and application of Density 
Functional techniques for use in molecular calculations. The complicated functional 
forms used in this theory preclude analytical integration, so that numerical methods 
must be used. Recently, Handy et a/.[l] proposed a radial integration based on the 
Euler-Maclaurin summation giving good results for a number of systems. However, 
it is worthwhile investigating other methods since for accurate calculations many 
quadrature points are required in general, even with the Euler-Maclaurin technique. 
In addition, on problem with the Euler-Maclaurin scheme is that there is little control 
over the location of additional quadrature points. For example, this scheme is very 
efficient for atoms, but in a molecular context additional points are required for an 
integration which is comparable to the atomic case. These additional points should 
ideally be placed in the valence region. However, with the Euler-Maclaurin scheme 
points must be placed throughout the entire region. This problem is bypassed with 
a multiregion integration. 

The advantage of the multiregion integration is the control of the placement 
of additional quadrature points in the molecular situation. In general, we determine 
an accurate quadrature for the constituent atoms or small molecular fragments. In 
the molecular situation we then only need to add points to the valence region(s), and 
increase the accuracy of the angular quadrature, in order to be able to produce an 
accurate molecular quadrature. One disadvantage of the multiregion integration is 
a possible increase in errors due to the accumulation of errors from the individual 
regions. 

Here we investigate the efficiency of a multi-region integration similar to that 
previously proposed by Te Velde and Baerends [2] and also Thakar [3]. One difference 
in the current work is the use of a change of variables due to McLean and Yoshimine [4] 
which allows for efficient placement of the quadrature points by mapping the region 
[a, m, 6] onto the region [—1, 0, 1], where a and b are the endpoints of the region, and 
m is any point between a and b. This mapping takes the form 

x = Cl + c 2 (i + /3)/(i - pt). 

The constants cj, c 2 and /? are completely determined by the mapping a — ► — 1, 
m — ► 0, and 6 — * 1. This mapping may be used for all ranges except the doubly 
infinite range (-oo, oo), providing appropriate limiting forms are used for special cases. 

The regions used here are based on the atomic charge density, in a spirit similar 
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to that used by Thakar [3]. However, rather than using a general rule based on the 
atomic number by which to determine the regions, we use the atomic charge density 
directly. The function r 2 * p(r) is generated numerically, and analysed for maxima, 
minima and points of inflection. In general, the regions we define are based on these 
points, in particular the points a and b map onto two minima in the charge density, 
and the point m maps onto the maximum in between these. These three points define 
a shell, and the mapping ensures an optimal distribution of quadrature points in the 
region. Generally we have one region for each shell, with the valence shell being split 
into an inner and outer parts at the point of inflection of the valence shell (see also 
Thakar), and the innermost (core) region into two at the first point of inflection. 
Although this scheme may seem somewhat complicated, in fact it is straightforward 
to generate the required regions and tabulate them for any atom. 

For the inner regions we use a standard Gauss-Legendre quadrature, whereas 
for the outer region we use a shifted Gauss-Laguerre scheme based on the weight 
function e"" r , where a is related to the highest occupied MO eigenvalue (or the IP) 
by a = 2(v/~ 2e t ) (see also Thakar [3]). 

The other details of the integration are very similar to that proposed by 
Handy et al. [1]. We use a product scheme for the angular integration identical 
to that of Handy et al., including the use of the angular crowding parameter Kg. The 
single center integration scheme of Becke [5] is used, but we use the cutoff function 
due to Handy et al., with standard Bragg-Slater radii. 

We consider three systems here - Ne, Cu and CO. We compare the Euler- 
Maclaurin to the current multiregion integration for all four systems. We consider 
integration of the total charge density only here, at the SCF level of approxima- 
tion, since Handy et al. [1] showed that this gives a good indication of the overall 
performance of an integration scheme. For Ne, we use the Dunning [6] [5s4p] con- 
tracted Gaussian basis set. For carbon and oxygen we use the correlation consistent 
polarized valence triple zeta (cc-pVTZ) basis sets of Dunning. For Cu we use a 
(20s 15p lOd 6/ 4p)/[6 + Is 5 + lp id 2/ lgr] basis set derived from the large primi- 
tive set of Partridge [7], contracted by Bauschlicher using the atomic natural orbital 
(ANO) approach. The C-0 bond length is 2.2 a.u. in all cases. 

The calculations were carried out using the Seward integral program and the 
Sweden SCF program on the Cray YMP-C90 at NASA Ames research center. 
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2 Results and Discussion 


For Ne atom, the results for the Euler- Maclaurin integration are given in Table 1. 
We give the number of points used in the radial integration, the number of points 
including those discarded due to the radial cutoff (in parentheses), the radial factor 
m T (Handy et al. [1], Eqn. 6), and the total error in the integrated charge density. 
Our results for this are very similar to Handy et al., as expected, with a very accurate 
radial integration attained with 72 radial points. In Table 2, we give the results from 
the multiregion integration, with the number of points in each region given. In region 
5 we also indicate the number of points including those discarded due to the radial 
cutoff, in parentheses. Overall, more radial points are required for the multi-region 
integration, although the integration is still accurate. For example, if one used a 
value for m r in the Euler-Maclaurin integration which was not optimal, then the 
results could be worse than for the multi-region integration. 

In Table 3, we give the results for the integration of Cu atom in a large ANO 
basis set, for three different numbers of radial points. One trend to be noted is that a 
higher m r value is needed for Cu than Ne as more radial points are added, and more 
radial points are needed to achieve a similar absolute accuracy. As noted by Handy et 
al. this is because the numerical integration only gives a certain relative error rather 
than absolute error in the integration. We note that for Cu atom (4s 1 ) it is necessary 
to integrate out a long way (22 a.u.) due to the very diffuse nature of the 4s orbital, 
and that sometimes a higher value of m r can be necessary than that recommended 
than Handy et al. to obtain optimum results. 

The multi-region results of Tables 4 and 5 compare quite well with the Euler- 
Maclaurin results, although again it seems that more points are necessary for the 
multi-region scheme, depending on the m T value used for in the Euler-Maclaurin 
scheme. We note that the inner regions for the multi-region scheme are very compact, 
due to the higher atomic number in this case, and that quite a few radial points are 
necessary to describe the density accurately in this compact region. However, this 
does not necessarily translate to a lot of points in a molecular calculation, since the 
number of angular points needed in the inner regions is much smaller than in the 
valence regions. 

Finally we consider the CO molecule (Tables 6-8). We use the same num- 
ber of angular points in for the Euler-Maclaurin and the multi-region integrations 
(n<?=42, n^=8 4), and the same Bragg-Slater radii in each case (rc=l. 32281 a.u., 
ro=l. 13383 a.u.). Overall, it seems that it is possible for the Euler-Maclaurin scheme 


4 



to outperform the multi-radial scheme if the right values of m r and m M (the angular 
factor [1]) are chosen. However, these differ significantly from the “standard” values 
recommended by Handy et al., m r = 2 or 3 and m M = 10 or 11. For these values, we 
see that the multi-region scheme is either equivalent to or better than the Euler- 
Maclaurin scheme in efficiency. Another interesting point from Tables 6 and 8 is that 
it is not just the total number of radial points which is important, but also the spread 
of these points along r, as this affects the total number of points through the angular 
factors. For example, the spread of radial points in the Euler-Maclaurin scheme is 
approximately linear on a logarithmic scale, apart from the very short and long range 
regions. Changing the factor m r changes the slope of this logarithmic plot, so that 
for higher values of m T there are more points in the core region and more points in 
the long range regions, leading to fewer points overall, since there are fewer angular 
points in the core region and the long range points axe discarded due to the radial 
cutoff. Thus for the Euler-Maclaurin scheme, very different total numbers of points 
are realized for m T = 3 versus m r = 4, even though the nominal number of radial points 
is the same in each case. 

For the multi-region scheme, two sets of results are presented for CO (Table 8), 
differing in the number of points in the valence and outer regions, and with the total 
number of points being very comparable to the best results found for the Euler- 
Maclaurin scheme. One advantage of the multi-region scheme is that we are able to 
take a set of integration parameters from a similar atom (for example, those for Ne in 
this case), and then place more radial points in the valence and outer valence regions 
in order to attain higher accuracy. This can be seen to be an effective way to add 
points for CO. 

Overall, it seems that neither integration scheme is clearly superior in the 
molecular situation, based on the current results. More and larger systems need to 
be studied in order to establish whether the multi-region integration scheme proposed 
here is significantly better to the widely used Euler-Maclaurin scheme of Handy et 
«/■[!] 
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Table 3: Euler- Maclaurin integration of Cu atom (r max =22.0 


n T 

m r 

Ap(r) 

47(64) 

2 

1 

O 

r*^ 

X 

© 

© 

43(64) 

3 

6.4xl0 -8 

40(64) 

4 

7.6xl0" 6 

71(96) 

2 

5.2xl0" 9 

64(96) 

3 

1.4xlO -10 

60(96) 

4 

4.4xl0 _1 ° 

95(128) 

2 

6.5xl0 -9 

86(128) 

3 

4.9xl0 -12 

80(128) 

4 

l.lxlO -13 

77(128) 

5 

l.lxlO -10 


Table 4: Regions for Cu atom (a.u.) 



a 

m 

b 

1 

0.000 

0.005 

0.010 

2 

0.010 

0.038 

0.077 

3 

0.077 

0.170 

0.380 

4 

0.380 

0.630 

0.940 

5 

0.940 

2.500 

6.000 

6 

6.000 

— 

22.000 
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Table 5: Multiregion integration of Cu atom (r max =22.0 a.u.) 


1 

2 

n T (region) 

3 4 5 

6 

Total n T 

Ap(r) 

15 

15 

15 

10 

20 

15(9) 

84(90) 

2.2xl0 -9 

20 

15 

15 

10 

20 

11(20) 

91(100) 

6.7xl0 -11 

25 

20 

20 

15 

25 

11(20) 

117(130) 

1.4xl0 -12 
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Table 6: Euler-Maclaurin integration of CO molecule (n$= 42, 

radbs=l. 32281, 1-13383 a.u.n^=84, r mar =15.0 a.u.). Note: same value of ra M and m T 
used for C and 0. 


n r 

Ngrid 

m r 


Ap(r) 

132,134 (192,192) 

60902,62666 

3 

10 

1.3xl0 -9 

132,134 (192,192) 

60902,62666 

3 

11 

3.0xl0 -1 ° 

132,134 (192,192) 

60902,62666 

3 

12 

2.8xl0 -11 

132,134 (192,192) 

60902,62666 

3 

13 

9.8xl0 -12 

132,134 (192,192) 

60902,62666 

3 

14 

2.5xlO -11 

148,150 (192,192) 

87162,88926 

2 

10 

2.1xlO -10 

148,150 (192,192) 

87162,88926 

2 

11 

8.7xl0 -u 

148,150 (192,192) 

87162,88926 

2 

12 

3.9xl0 -11 

148,150 (192,192) 

87162,88926 

2 

13 

2.2xl0 -11 

148,150 (192,192) 

87162,88926 

2 

14 

1.6xl0 -11 

148,150 (192,192) 

87162,88926 

2 

15 

1.7xl0 -11 

124,125 (192,192) 

47228,48110 

4 

11 

8.7xl0 -9 

124,125 (192,192) 

47228,48110 

4 

12 

5.4xl0 -9 

124,125 (192,192) 

47228,48110 

4 

13 

1.8xl0 -9 

124,125 (192,192) 

47228,48110 

4 

14 

5.2xl0 _1 ° 

124,125 (192,192) 

47228,48110 

4 

15 

1.5xl0 -9 
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Table 7: Regions for CO molecule (a.u.) 



a 

m 

b 

1 

0.00 

c 

0.025 

0.05 

2 

0.05 

0.17 

0.63 

3 

0.63 

1.27 

1.90 

4 

1.90 

4.0 

6.0 

5 

6.0 

— 

15.0 

1 

0.0 

0 

0.02 

0.04 

2 

0.04 

0.12 

0.41 

3 

0.41 

0.83 

1.37 

4 

1.37 

4.0 

6.0 

5 

6.0 

— 

15.0 






Table 8: Multiregion integration of CO molecule (r max 15.0 a.u.) 




n r ( region) 


Total n r 



A/>(r) 

1 

2 

3 

4 

5 





mm 

MM 

24 

24 

7(16) 

90(99) 

58078,56944 

12 

3.0xl0~ 9 

K9 

m 

24 

24 

7(16) 

90(99) 

58078,56944 

13 

l.OxlO -10 

15 

20 

24 

24 

7(16) 

90(99) 

58078,56944 

14 

1.4xl0 -9 

15 

20 

24 

24 

9(24) 

92(107) 

59842,58708 

13 

8.3xl0 _ " 

15 

20 

24 

24 

11(32) 

94(115) 

61606,60472 

13 

8.3xl0~" 

15 

20 

24 

32 

9(24) 

100(115) 

66898,65764 

13 

5.7x10-" 

15 

20 

24 

40 

9(24) 

108(123) 

73954,72820 

13 

7 

o 

r— H 

X 

o 

15 

20 

24 

40 

9(24) 

108(123) 

73954,72820 

12 

1.4xl0" u 

15 

20 

24 

40 

9(24) 

108(123) 

73954,72820 

14 

1.2x10-" 
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